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The extension of Boltzmann-Gibbs thermostatistics, proposed by Tsallis, introduces an additional 
parameter q to the inverse temperature j3. Here, we show that a previously introduced generalized 
Metropolis dynamics to evolve spin models is not local and does not obey the detailed energy 
balance. In this dynamics, locality is only retrieved for q — 1, which corresponds to the standard 
Metropolis algorithm. Non-locality implies in very time consuming computer calculations, since the 
energy of the whole system must be reevaluated, when a single spin is flipped. To circumvent this 
costly calculation, we propose a generalized master equation, which gives rise to a local generalized 
Metropolis dynamics that obeys the detailed energy balance. To compare the different critical values 
obtained with other generalized dynamics, we perform Monte Carlo simulations in equilibrium for 
Ising model. By using the short time non-equilibrium numerical simulations, we also calculate for 
this model: the critical temperature, the static and dynamical critical exponents as function of 
q. Even for q 7^ 1, we show that suitable time evolving power laws can be found for each initial 
condition. Our numerical experiments corroborate the literature results, when we use non-local 
dynamics, showing that short time parameter determination works also in this case. However, the 
dynamics governed by the new master equation leads to different results for critical temperatures 
and also the critical exponents affecting universality classes. We further propose a simple algorithm 
to optimize modeling the time evolution with a power law considering in a log-log plot two successive 
refinements. 



PACS numbers: 05.10.Ln, 05.70.Ln, 02.70.Uu 

I. INTRODUCTION 

The study of the critical properties of magnetic sys- 
tems plays an important role in statistical mechanics and 
as a consequence also in thermodynamics. For equilib- 
rium, the extensitivity of the entropy is a question of 
principle for most physicists. Nevertheless, an important 
issue may be raised. While many physicists believe that 
statistical mechanics generalizations with an extra pa- 
rameter q |l[ are suitable to study the optimization com- 
binatorial process as for example the simulated anneal- 
ing (see e.g. or areas sucn as econophysics 0,0, 
population dynamics and growth models [6|-|9(, Bibliom- 
etry [Io| and others. 
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In this paper, we generate the critical dynamics of 
Ising systems using a new master equation. This master 
equation leads to a generalized Metropolis prescription, 
which depends only on the spin interaction energy vari- 
ations with respect to its neighborhood. Furthermore, 
it satisfies the detailed energy balance condition and it 
converges asymptoticall y to the generalized Boltzmann- 
Gibbs weights. In Refs. |TD.[l2j generalized prescriptions 
have been treated as local. Here, we demonstrate that 
they are instead non-local. However, a non-local pre- 
scription such as the one of Ref. [l3j is numerically more 
expensive and destroys the phase transition. Another 
possibility is to recover locality. Using a special defor- 
mation of the master equation, we show how to recover 
locality for a generalized prescription and additionally re- 
covering the detailed energy balance in equilibrium spin 
systems, maintaining the system phase transition. 

To apply our Metropolis prescription, we have sim- 
ulated a two dimensional Ising system in two different 
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ways: using equilibrium Monte Carlo (MC) simulations 
we estimate critical temperatures for different q-values 
and performing time-dependent simulations. In the sec- 
ond part, we also calculate the critical exponents set cor- 
responding to each critical temperature. Finally, we have 
developed an alternative methodology to refine the de- 
termination of the critical temperature. Our approach 
is based on the optimization of the magnetization power 
laws in log scale via of maximization of determination 
coefficient (r) of the linear fits. 

Our presentation is organized as follows. In Sec. [Til 
we briefly review the results of the critical dynamics for 
spins systems. In this review, we calculate the critical 
exponents for the several spin phases, that emerge from 
different initial conditions. In Sec. IIII1 we propose a new 
master equation that leads to a Metropolis algorithm, 
which preserves locality and detailed energy balance, also 
for g^l. In Sec. IIV1 we simulate an equilibrium Ising 
spin system in a square lattice and show the differences 
between the results of our approach and of Refs. (lHll^. 
Next, we evolve a Ising spin system in a square lattice, 
from ordered and disordered initial conditions in the con- 
text of time dependent simulations. From such non equi- 
librium Monte Carlo simulations, also called short time 
simulations, we are able to calculate the dynamic and 
static critical exponents ones. Finally, the conclusions 
are presented in Sec. [V] 



II. CRITICAL DYNAMICS OF SPIN SYSTEMS 
AND TIME DEPENDENT SIMULATIONS 

Here, we briefly review finite size scaling in the dy- 
namics relaxation of spin systems. We present our al- 
ternative deduction of the some expected power laws in 
the short time dynamics context. Readers, which want 
a more complete review about this topic, may want to 
read 

This topic is based on time dependent simulations, and 
it constitutes an important issue in the context of phase 
transitions and critical phenomena. Such methods can 
be applied not only to estimate the critical parameters 
in spin systems, but also to calculate the critical expo- 
nents (static and dynamic ones) through different scaling 
relations by setting different initial conditions. 

The study of the statistical systems dynamical crit- 
ical properties has become simpler in nonequilibrium 
physics after the seminal ideas of Janssen, Schaub and 
Schmittmann [l5j and Huse [l6j . quenching systems from 
high temperatures to the critical one, they have shown 
universality and scaling behavior to appear already in 
the early stages of time evolution, via renormalization 
group techniques and numerical calculations respectively. 
Hence, using short time dynamics, one can often circum- 
vent the well-known problem of the critical slowing down 
that plagues investigations of the long-time regime. 

The dynamic scaling relation obtained by Janssen et 
al. for the magnetization fc-th moment, extended to finite 



size systems, is written as 

(M k )(t,T,L,m ) = ^/"(M^r^^r'L.^mo), 

(1) 

where the arguments are: the time t; the reduced temper- 
ature t = (T — T c )/T c , with T c being the critical one, the 
lattice linear size L and initial magnetization too. Here, 
the operator (. . .) denotes averages over different configu- 
rations due to different possible time evolution from each 
initial configuration compatible with a given ttiq. On the 
equation right-hand-side, one has: an arbitrary spatial 
rescaling factor 6; an anomalous dimension xq related to 
mo- The exponents /? and v are the equilibrium critical 
exponents associated with the order parameter and the 
correlation length, respectively. The exponent z is the 
dynamic one, which characterizes the time correlations in 
equilibrium. After the scaling 6 _1 L = 1 and at the criti- 
cal temperature T = T c , the first [k = 1) magnetization 
moment is: (M)(t,L,m ) = L~^ u (M)(L- Z t, L x °m ). 

Denoting u = tL~ z and w — L x °mo, one has: 
(M)(u,w) = (M)(L~ z t,L x °m ). The derivative with 
respect to L is: d L {M) = {-p/v)L-Pf u - x {M)(u,w) + 
L~P/ v [d u (M)dLU + d w (M)dLw], where we have explic- 
itly: 8lu — —ztL^ 2 ^ 1 and &lw = xomoL x °~ 1 . In the 
limit L — > co, 8l(M) — > 0, one has: xowd w (M) — 
zud u (M) — j3/v(M) = 0. The separability of the vari- 
ables u and w in (M)(u,w) = Mi(u)M 2 (w) leads to 
X0WM2/M2 = P/v + zuM[/M2, where the prime means 
the derivative with respect to the argument. Since 
this equation left-hand-side depends only on w and 
the right-hand-side depends only on u, they must be 
equal to a constant c. Thus, Mi(u) = u^/*)-^"*) 
and M^iw) — w c ^ x °, resulting in (M) (u,w) — 
m c J x a ■£ j p/uj.(c-p/v)/z _ Returning to the original variables, 

one has: (M)(t, L, m ) = rn^t^l^- lz . 

On one hand, choosing c = xq and calculating 9 = 
(xq — P/v)/z, at criticality (r = 0), we obtain (M) mo ~ 
mot corresponding to a regime under small initial mag- 
netization. This can be observed by a finite time scaling 
b = i 1 / 2 in equation [TJ at critical temperature (t = 0) 
which leads to (M) (t,m ) = t~ f3 ^ v ^{M){l,t x °/ z m ). 
Defining x = t x °' z mo, an expansion of the averaged 
magnetization around x — results in: (M)(l,x) = 
(M)(1,0) + d x (M)\ x=Q x + 0{x 2 ). By construction 
(M)(1,0) = 0, since u = t x °/ z m a < 1 and d x {M)\ x=0 
is a constant. So, by discarding the quadratic terms we 
obtain the expected power law behavior (M) mo ~ mot 6 . 
This anomalous behavior of initial magnetization is valid 
only for a characteristic time scale i max ~ Tn z ^ x ° . 

On the other hand, the choice c = corresponds to 
a case where the system does not depend on the initial 
trace of the system; and Too = 1 leads to simple power 
law: 

(M) roo=1 ~ t-VW (2) 

that similarly corresponds to decay of magnetization for 
t > tmax of a system that previously evolved from a ini- 
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tial small magnetization (mo), and had its magnetization 
increased up to a magnetization peak. 

For m = 0, it is not difficult to show that the magne- 
tization second moment is 

(M 2 ) n ~ t«-WMI» ( 3 ) 

\ /mo— v ' 

where d is the system dimension. 

Using Monte Carlo simulations, many authors have ob- 
tained the dynamic exponents 9 and z as well as the 
static ones /3 and v, and other specific exponents for 
many different models and situations: Baxter- Wu (l7j . 
2, 3 and 4-state Potts [l8], fl9| . Ising with multispin 
interactions [2(|, models with no defined Hamiltonian 
(celular automata and contact process) |2lT[23l|. models 
with tricritical point (24|, Heisenberg j25||, protein fold- 
ing [H, , propagation of damages in Ising models [28[ . 

The sequence to determine the static exponents from 
short time dynamics is: to determine z first, perform- 
ing Monte Carlo simulations that mixes initial condi- 
tions [18j . and consider the power law for the cumulant 

( M2 ) n ,/ 

w)= w 
{ M / mo =i 

Once z is calculated, the exponent 77 = 2/3/ f is cal- 
culated according to r\ — 2(/3/pz) ■ % where (ft/vz) was 
estimated via magnetization decay and z from cumulant 
F 2 . 

However, prior to obtaining the critical exponents, we 
also perform time dependent MC simulations in order 
to refine the critical temperatures. These are based on 
power laws obtained by finite size scaling analysis of 
the magnetization decay from an initially ordered state 
(Eq. This choice demands a number of runs smaller 
than other power laws in non-equilibrium, and so we pro- 
pose an simple algorithm that spans different critical val- 
ues to find the best determination coefficient in linear fit 
ln(M) versus \nt . This procedure is explored in Sec. lIVl 
and is used later to calculate the critical temperatures for 
Ising models with different values of the non-extensivity 
parameter q in our new Metropolis prescription. 



A. Standard master equation and Metropolis 
algorithm 

In general, spin systems non-equilibrium dynamics 
are described by the time evolution of the probability 
P(E,t) that, at instant t, the system has an energy E. 
This probability is obtained from the master equation: 
dP{E^ a \t)/dt = £{Maf° -> a\ a) ]P[E( b \t]-w[a\ a) -> 

af']P[E^ a \ t]}, where w[uf' — > a^} is the transition 

rate of the i-th spin from erf } to a[ a) . Here, E^ (E^) 
is the energy of the system before (after) the transition. 
As t — > 00, dP(E,t)/dt — is a necessary condition 
for equilibrium. A sufficient but not necessary condi- 
tion for equilibrium, known as detailed balance condi- 
tion, supposes a more restricted situation for ocurrence 
of dP(E,t)/dt = 0, i.e., w[af ] -> a\ a) ]P[E^]-w[(x\ a) -> 
a^']P[E^] = 0, meaning that each term in the summa- 
tion vanishes. In this case, P(E) = P(E, t — > 00) is 
the Boltzmann distribution: P(Ej) = e^ B >/Et e~ pEk , 
where the summation is over the different energy states 
and (3 = (fcsT)- 1 . 

Employing detailed balance requires to find simple 
prescriptions for spins systems dynamics, as for ex- 
ample the Metropolis prescription: w[<T( — > (r'f ] = 
min{l, exp[— (3(E^ — E^)}}. When applied to evolve 
spin systems, this simple dynamics reduce to calculate 
just local energy changes. For instance, the Ising model 
in two dimensions has an energy E^> — —Ja^' i Si x t i y +£ 
before the flip of spin <7< a) t , located at site indexed by i x 
and i y , where the local energy change is quantified by 

Si x ,i y = Gix + l.iy + X,i y + Qi^.iy—l + ^i^.iy + l 

and the non-local energy is £, which is obtained excluding 
the spin Ot S) i , from the calculation. After the spin flip, 

the energy is E^ = —Jo'i^S^^ + £ and the energy 
change of the system, due to the spin cr^i flip is simply: 

E^-E^ = -J\aQ iy -a^ iy \S ix , iy , (5) 
which does not depend on the energy of the other spins. 



III. GENERALIZED MASTER EQUATION 

In this section, we start recalling the way that the 
Metropolis algorithm is obtained from the master equa- 
tion for spin systems. We point out that the energy dif- 
ference by flipping an Ising spin is local, i.e. it depends 
only on the flipped spin. Next, we show a first attempt 
to generalize the Metropolis algorithm [TJ [l2| , accord- 
ing to the non-extensive thermostatistics, introduced by 
Tsallis 1.]. We show that this generalization does not 
preserve the spin flip locality. To recover this locality, we 
propose a new generalized master equation, which leads 
to a different generalization of the Metropolis algorithm. 



B. Generalized Metropolis algorithm 

The system equilibrium is described by the generalized 
Boltzmann-Gibbs distribution 



Pi- g (Ei) 



E?=i[ei- 9 (-/3'£i)]« 



(6) 



where f2 is the number of accessible states of the sys- 
tern and /3' = /3/ £{[ei-,(-«] 9 + (1 - g)/3(£)i-J, 



i=l 



where (E)\- q = £ i=1 EiPi- q (Ei). Here it is important 
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to mention that (fcs/?') -1 is a scale temperature that can 
be used to interpret experimental and computational ex- 
periments. There is a heated ongoing discussion whether 
it is the physical temperature or not. 
The function 



«0) = 



_ / (1 + aar) 



l/a 







for ax > - 
otherwise 



(7) 



is the generalized exponential [30, [3l[. For a — > 0, one 
retrieves the standard exponential function en(x) = e x . 
It is this singularity at ax > — 1 that brings up inter- 
esting effects such the survival/extinction transitions in 
one-species population dynamical models [<S]. The in- 
verse of the generalized exponential function is the gen- 
eralized logarithmic function \n a (x) = {x a — l)/ot, which 
for a — > leads to the standard logarithm function 
lno(x) = ln(a;). Notice that the inequality ax > — 1, 
for fixed x produces a limiting value for a. This gen- 
eralized logarithmic function has been introduced first 
in the context on non-extensive thermostatistics 0, |30] 
and has a clear geometrical interpretation ss the area 
between 1 and x underneath the non-symmetric hyper- 
bole \jt x ~ a [3l[. It is interesting to notice, that in 1984 
Cressie and Read [32[ proposed an entropy that would 
lead to a generalization of the logarithm function given 
by : ki a (x)/(a + 1). In this case, we would gain the lim- 
iting value in a but lose its geometrical interpretation. 

To recover the additive property of the argument, 
when multiplying two generalized exponential functions: 
e a (a)e a (b) — e a (a ® a b) [e a (a) / e a (b) — e a (a Q a b)] and 
e a (a) (g) Q e a (b) — e a (a + b) [e a (aj Q) a e n (b) = e n (a — b)] 
consider the following algebraic operators [33|, [34| : 



a ©c b = a + b + aab 
a — b 

a Q a b = 



1 + ab 

a® a b = (a a +b a -l) 1/a 
aQ a b= (a a -b a + l) 1/a 



(8) 
(9) 
(10) 
(11) 



Observe that, if a Q a 6 = 0, then a = b and if a <Q a b — 
c <g) a d, then a Q c — d Q a b. 

However, in equilibrium, the Ising model prescribes an 
adapted Metropolis dynamics that considers a general- 
ized version of exponential function [IJ [l2| : 



r 0>) i a • 



(o)l _ Pi- q [E^] _ <e^ q [-P'E^] 



Pi- q [E(»)} 



ei- q [ 



_ . (12) 

From the generalization of the exponential function 
in the Boltzmann-Gibbs weight, the transition rate of 
Eq.HUcan be used to determine the system evolution, as 

I 



the Metropolis algorithm. Nevertheless, we stress that in 
such a choice, the dynamics is not local. Because gener- 
alized exponential functions are non-additive, a spin flip 
introduces a change in the system energy that is spread 
all over the lattice. More precisely, consider the Ising 
model in a square lattice, one can show that: 

^M — -t-^ 6 -^" (13) 



where E^ - E^ is given by Eq. [51 which depends only 
the spins that directly interact with the flipped spin, vi- 
olating the detailed energy balance. 

In Refs [ll|, [T3| , the authors consider (with no expla- 
nations) the equality in Eq. Q31 instead of considering 
Eq. [13] Thus, the detailed energy balance is violated, 
since the system is updated following a local calculation 
of the generalized Metropolis algorithm of Eq. [T2] 

To correct this problem, one must update the spin sys- 
tem using the non- locality of Eq. 1121 which is numerically 
expensive, since the energy of the whole lattice must be 
recalculated due to a simple spin flip. The other alterna- 
tive is to require that the transition rate depends locally 
in the energy difference of a simple spin flip, which in 
turn leads us to a modified master equation. Since the 
former is very expensive numerically, we explore only the 
latter alternative which is numerically faster and is able 
to produce statistically significant results for fairly large 
spin systems. 



C. Recovering locality in the generalized 
Metropolis algorithm 

Based on the operators of Eq. [8] to Eq. [Til we propose 
the following generalized master equation: 



dPt-gjE^ 

dt 



E 



i/'IffW ^^ a) ]Q q/q P q [E^] Q d 



q/q 



w[al a) ^^ b) ]Q q/q P q [E^}. (15) 

where P q (E) is given by Eq. [5] Here, it is suitable to 
call q = 1 — q and write the generalized exponentials 
as a function of q. In equilibrium, dP\- q /dt = and a 
dynamics governed by Eq[Bl 

The detailed balance (a sufficient condition to equilib- 
rium) for the generalized master equation is 



h ®q,qW[4 a) 



— ¥ a 



(6)1 



P q [E^]Q q/q P q [E^}, 



(16) 
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which leads to a new generalized Metropolis algorithm: 



(<>) 



o) ) =min{l, [egi-P'EW)]* 2> 9/q [e 9 (-(3'E^)} "} = min {l, [e 9 ~(-/3'(£ (a) - E®)j\ "} 

J6) 



iin{l, [e$(/3' 



J 



r («) 



'} 



(17) 



and now the transition probability depends only on en- 
ergy between the read site and its neighbors, i.e., locality 
is retrieved. 



IV. GENERALIZED METROPOLIS 
ALGORITHM - NUMERICAL SIMULATION 
RESULTS 

We have performed Monte Carlo simulations of the 
square lattice Ising model in the context of generalized 
Boltzmann-Gibbs weights. These simulations are based 
on two approaches for Metropolis dynamics. The first 
one (Metropolis I) is described in Ref. [ll|, where the 
nonlocal transition rate of Eq. ([T2|) is used to update 
the spin system. In the second approach (Metropolis II) , 
the local transition rate of Eq. [17] is used. We separate 
our results in two different subsections: the equilibrium 
simulations and short time critical dynamics. 



A. Equilibrium 

In this part we analyze the magnetization (to) , where 
(•) denotes averages under Monte Carlo (MC) steps. We 
perfom MC simulations for q = 0.6, q = 0.8 and q = 1.0. 
In the simulations, we have used L m i n = 2 4 = 16 up 
to L max = 2 9 = 512, with periodic boundary condi- 
tions and a random initial configuration of the spins with 
(too) = 0. Differently from reported in Ref. [Til where the 
results have been obtained after 10 7 MC steps per spin, 
we have used 6.1 3 MC steps per spin, an equilibrium sit- 
uation consistent with the one reported by Newman and 
Barkema [Hj]. This results in 1.5 ■ 10 6 - 1.5 • 10 9 MC steps 
for the whole lattice of 16 2 up to 512 2 spins. 

Fig. [T] shows the magnetization curves as function of 
critical temperature for different q-values. The critical 
temperature increases as q decreases. This behavior, us- 
ing our algorithm (Metropolis II) differs from the one ob- 
tained using the algorithm of Refs. [ll| and[l2| (Metropo- 
lis I). We stress that both algorithms agree for q = 1, 
the usual Boltzmann-Gibbs weights, converging to the 
theoretical value log(l + y/2)/2. In Table [fl we show 
the critical temperature and error obtained from the ex- 
trapolation L — > oo (see Fig. [5]) using both algorithms. 
These results suggest a thorough difference among the 
processes and critical values found between two the dy- 
namics Metropolis I and II. In Fig. Q] the curves show 
phase transitions for critical values upper to log(l+v / 2)/2 
as q < 1. This differs from previous studies, which are 
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FIG. 1. System magnetization versus temperature for q — 

I. 0, 0.8 and 0.6. Using the dynamics based on Metropolis 

II, we observe phase transitions for critical values upper to 
log(l + v / 2)/2 as q < 1 differently from previous studies, which 
are based on Metropolis I. 



based on prescription Metropolis I. 

Fig. Q] shows that, differently from the q — 0.8 and 
q = 1.0 cases, for q = 0.6 the discontinuity in the mag- 
netization curve does not depend on system size L. In 
fact, in this case, the critical temperature T c does not 
depend on L. This effect occurs due to the cutoff of the 
escort probability distribution as reported for Metropolis 
I [lH for q < 0.5. For Metropolis II, Fig. H depicts that 
T c remains constant for all values of L~ l , for q = 0.6. 
For both cases, q = 1.0 (obviously) and q = 0.8, we have 
verified that v w 1 and j3 « 0.125, obtained from the col- 
lapse of the curves (M) L^l u versus (T-T C )L^ V . This 
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data collapse permits the extrapolation of k B T c / J ver- 
sus L" 1 , since v « 1 for both cases according to Fig. 
In following, we show using non-equilibrium simulations 
that 2/3 /V « 0.25, for g = 1 and g = 0.8, validating the 
data collapse results (see table [Vj) . 




0.00 0.02 0.04 0.06 0.08 0.10 
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FIG. 2. Extrapolation (L — > oo) of critical temperatures for 
different q— values: 0.6, 0.8 and 1.0 for the 2d Ising model. 

Another important question to be formulate is: Can we 
corroborate the same behavior in non-equilibrium simu- 
lations? Next section, we show results from MC sim- 
ulations in non-equilibrium regime under the two dy- 
namics (Metropolis I and II). We also analyze the criti- 
cal exponents (dynamic and static) as a function q from 
short-time dynamics. We show that short time dynam- 
ics corroborate the behavior predicted by two dynam- 
ics suggesting that Metropolis II indeed presents an in- 
crease of critical value as q- value increases different from 
Metropolis I. Our results suggest that these technics 
based on timc-dcpcndcnt simulations can be extended 
also for q ^ 1, in short range spin models. 



B. Short time 

Here we address time dependent MC simulations in 
the context of so called short time dynamics. First, to 
test our methodology, we show that critical values ob- 
tained from non-equilibrium simulations using Metropo- 
lis I must corroborate the critical values obtained in 
Ref. [Ill ], where MC simulations at equilibrium have 
been employed. We have checked it. Nevertheless, as 
in the equilibrium numerical simulations, we show that 
Metropolis II leads to different values from Metropolis I 
method. 

Our algorithm to estimate the critical temperature is 
divided in two stages. In the first stage, a coarse grained 
calculation is performed to estimate the critical tempera- 
ture T c (q), for different q values. In the second stage, one 
uses the estimated critical temperature obtained in the 



first stage to run a non-equilibrium Monte Carlo simula- 
tion. We denote the second state as fine scale stage. In 
this stage, one determines the dynamical critical expo- 
nent from the short time behavior of the spin system, as 
described in Sec. [TIJ Since, even using non-extensive ther- 
mostatistics, the magnetization must behave as a power 
law (M) ~ t~Pl vz , we conjecture that changing T c (q) 



from T c (min) (<?) up to T^'(q), the best T c (q) is the one 
that leads to the best linear behavior of ln(M) versus 
Int. We have considered n s = 500 realizations, with ini- 
tial magnetization mo = 1. 

From the theoretical critical temperature (f3 c = 
J/k B T c = log(l + V2)/2), one allows the temperature 
to vary in the range from k B T c / J — 1 up to k B T c / J + 1, 
setting k B T/J = [2 - log(l + V2)][log(l + >/§)] + J • A, 
where A = 0.1 and j = 0,1,..., 20. This is the coarse 
grained stage. For each temperature, a linear fit is per- 
formed and one calculates the determination coefficient 
of fit as: 



-,(max) 



Nmc 

£ (In (A/) -a~b\ntf 
t=i 

Nmc 



(18) 



E (hr(M)-ln(M)(i))2 
t=i 



and ln(Af) = (1/N MC ) E t =T MM) (*)> where N mc is 
the number of Monte Carlo sweeps. In our experiments, 
we have used Nmc = 300 MC steps. Here, r = 1 means 
an exact fit, so that the closer r is from the unity, the bet- 
ter. Here, a and b are the linear coefficient and the slope 
in the linear fit ln(Af) versus hut, respectively. From b, 
one estimates the exponent —/3i//z. 

In the fine scale stage, we refine the critical temper- 
ature k B Tc 1 \q)/J obtained in the first stage. We use 
the algorithm considering A — 0.01, with j = 0, 1, . . . , 20 
considering k B T^\q)/J = k B T i c 1 \q)/J - 0.1 + j ■ A, 
now to find the best critical temperature in the range 
k B TP(q)/J - 0.1 to k B TP(q)/J + 0.1 with precision 
A = 0.01. 

A natural validation for our algorithm is to reproduce 
the results obtained in Ref. [ill in equilibrium, using the 
Metropolis I approach, for a specific q value, considering 
our MC non-equilibrium simulations. For instance, for 
q = 0.70, one has at equilibrium k B T c /J = 1.891(7), 
in Ref. [ll(. After two stages, our algorithm produces 

k B Tc / J = 1.889, validating our numerical code. 

Next, we use the algorithm with the following values: 
q = 0.70, 0.75, 0.80, 0.85, 0.90, 0.95 and 1.00, in the equi- 
librium situation. In Table UH we show our results for the 
first stage (coarse grained) using Metropolis I prescrip- 
tion. The values of the determination coefficient a of the 
linear fit ln(Af) versus hit are presented for different q 
values. The highest values (in bold) correspond to best 
critical temperature found in the first stage. For example: 
for q = 0.75, we have that best a value is 0.940558872, 
which corresponds to k B T^\q)/J = 1.86918531. 

In Table m the symbol "-" corresponds to situations 
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q ref. Metropolis I Metropolis II 
0.6 1.761(3) 3.201(1) 
0.8 1.891(7) 2.461(5) 

1.0 2.259(11) 2.262(9) 



TABLE I. Comparison between critical the critical temperature and error, for the 2d Ising model, obtained from extrapolation 
L — > oo (see Fig. [2| using the algorithm of Ref [Tl|] (Metropolis I) and ours (Metropolis II) . 



where the computation of slopes is not possible, due to 
large deviations in magnetization. 

After the refinement (second stage), the best values 
found for the critical temperatures using Metropolis I 
prescription, for different q— values are presented in the 
first line of Table EH In Fig. El for q = 0.70 and 0.85, we 
show the magnetization decays as the power law: M(t) ~ 
t~P/ vz , for the critical temperature estimated using: our 
(Metropolis II) and Metropolis I algorithms. Also, we 
show the plots considering MC simulations for T c + 6 
and T c - S, with S = 0.05. 



o 

0.6 




A T - 0.05 
Generalized Metropolis I 

t(MCSteps) 




A T - 0.05 
Generalized Metropolis I 



t(MCSteps) 



1 0.8; q = 0.70 



01 

C0.6 

Ml 

rtj 

0.4- 




A T - 0.05 
Generalized Metropolis I 



t (MCsteps) 




10 100 

t(MCsteps) 



FIG. 3. Decay of magnetization according to the power law: 
M(t) ~ t~Pl vz i n the critical temperature found by the con- 
sidered algorithms (circle red points), for q — 0.70 and 0.85. 
We also show the plots considering MC simulations for T c + 5 
and T c — 8. It was used 5 = 0.05. The upper (lower) plots 
correspond to Metropolis I (II) algorithm. 

We use the same procedure to find the critical temper- 
atures for prescription Metropolis II. We find very differ- 
ent results, when compared with that ones obtained with 
Metropolis I. Similarly to Table HH we show the results 
using the Metropolis II prescription in Table IIVI The 
values are smaller than the ones found with Metropolis 

I prescription. However, they match as q — > 1, which 
validates the numerical procedure. 

Similarly, the best results after the fine scale refinement 
(second stage) are shown in the first line of Table [Vj 
The magnetization decay obtained by the Metropolis 

II algorithm depicted in Fig. [31 After obtaining these 
estimates for the critical temperatures, we perform short 



time simulations to obtain the critical dynamic exponents 
z and 9 and the static one rj = using the power laws 

of Sec. [TTJ Here, we calculated 8 from time correlation 
C(t) = (M(i)Af(O)). Tome and de Oliveira [H showed 
that correlation behaves as C (t) ~ t e , where 9 is exactly 
the same exponent from initial slope of magnetization 
from lattices prepared with initial fixed magnetization 
nio- The advantage of this method is that we repeat 
N s runs, but the lattice does not require a fixed initial 
magnetization. It is enough to choose the spin with prob- 
ability 1/2. i.e., mo = in average. This method does 
not require the extrapolation too — > 0. 

Figs. H] and [5] depict plots of time evolving of F 2 of Eq.@] 
and C(t) as a function of t for the different Metropolis 
algorithm. 

To obtain the exponents consider the following steps. 
Firstly, in simulations that start from the ordered state 
too = 1 and L — 512, calculate the slope j3/i/z of the 
linear fit of \n(M(t)) as a function of Int. The error 
bars are obtained, via running simulations for Nbi„ = 5, 
calculating (M(t)) that each for seed, with N run = 400 
runs. 



q = 0.70 



Generalized Metropolis I 



q = 0.85 



Generalized Metropolis I 



q = 0.70 




t(MCSteps) 



Generalized Metropolis II 



q = 0.85 




Generalized Metropolis I 



10 100 

t(MCsteps) 



10 100 

t(MCsteps) 



FIG. 4. Dynamic cumulant Fa(t) versus t in log scale. The 
slope gives d/z which supplies the z- value. Both prescriptions 
(Metropolis I e II) are studied. 

Once we have calculated fi/vz, we estimate z tak- 
ing the slope in log-log plot I11F2 versus Int. We used 
N s = 3000 different runs starting from random spins con- 
figurations with too = for time series (M(t) 2 ) x t and 
the same number of runs for time series (M(t)) x t start- 
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k B T^/J 


g = 0.70 


q = 0.75 


q = 0.80 


q = 0.85 


q = 0.90 


q = 0.95 


q = 1.00 


T* — 0.6 


644974839 


544533372 

V / . ' ) i t')')' J') 1 


433158672 


376097284 


0.326410245 


0.282451042 


263459517 


T* — 0.5 


998967222 


872350719 


65553425 


492266999 


392016875 


33691 5743 


306414971 


T* — 0.4 


858060019 


0.940558872 


0.979041762 


731374836 


500535816 


396550709 


341806493 


T* - 0.3 


0.822853648 


0.82063843 


0.90326648 


0.999101612 


0.773207193 


0.535343676 


0.409044416 


T* - 0.2 






0.833548876 


0.885834994 


0.998950355 


0.788883669 


0.547803953 


T* - 0.1 








0.836458324 


0.882565862 


0.999817616 


0.776353951 


T* = log(l + v/2)/2 










0.817075651 


0.897612219 


0.997114577 


T* +0.1 












0.82435859 


0.916225167 



TABLE II. Coarse grained Stage for Metropolis I. The values of determination coefficient a of the linear fit ln(M) versus Int 
for different q values. The highest values are in bold and correspond to best critical temperature found at first stage (coarse 
grained). For example: for q=0.75, the best r is 0.940558872, which corresponds to kaTc^ / J=l. 86918531. 



q 


0.70 


0.75 


0.80 


0.85 


0.90 


0.95 


1.00 




1.77(1) 


1.82(1) 


1.89(1) 


1.97(1) 


2.07(1) 


2.17(1) 


2.27(1) 


p/uz 


0.060(4) 


0.062(7) 


0.078(5) 


0.082(2) 


0.100(5) 


0.094(4) 


0.057(3) 


z 


2.13(4) 


2.15(5) 


2.12(4) 


2.09(3) 


2.10(3) 


2.11(6) 


2.15(3) 


e 


0.18(4) 


0.14(4) 


0.22(7) 


0.17(3) 


0.04(6) 


0.17(3) 


0.19(4) 


V 


0.25(2) 


0.27(3) 


0.33(2) 


0.34(1) 


0.42(2) 


0.40(2) 


0.25(1) 


r 


0.998568758 


0.998915473 


0.999342437 


0.999458152 


0.999589675 


0.999718708 


0.999206853 



TABLE III. Critical temperature and exponents obtained for different q values for prescription Metropolis I. The exponents 
where obtained performing simulations for the estimated critical temperatures and were based on power laws previously de- 
scribed in short time regime. The last line we show the r value for the best fits in the second stage (fine scale) 





q = 0.70 


q = 0.75 


q = 0.80 


q = 0.85 


q = 0.90 


q = 0.95 


q = 1.00 


T* - 0.1 




0.369354816 


0.393077206 


0.473151224 


0.555233203 


0.664145691 


0.754299904 


T* = log(l + \/2)/2 




0.439789081 


0.50595381 


0.658725599 


0.829581264 


0.952694449 


0.997206853 


T* +0.1 




0.579484235 


0.756653731 


0.959530136 


0.995694839 


0.951627799 


0.91284158 


T* +0.2 


0.601895662 


0.836896384 


0.999315534 


0.932708995 


0.869271327 


0.848547425 


0.838519074 


T* +0.3 


0.844382176 


0.989767198 


0.875131078 


0.833730495 


0.831169795 


0.799709762 




T* +0.4 


0.989716381 


0.847004746 


0.812106454 










T* +0.5 


0.828738110 


0.787203767 












T* + 0.6 


0.842827863 















TABLE IV. Coarse grained Stage for Metropolis II - The values of determination coeficient r of the linear fit ln(M) versus lnt, 
for different q values. As in Table [TO the highest values are in bold correspond to best critical temperature found at first stage 
(coarse grained). 



q 


0.70 


0.75 


0.80 


0.85 


0.90 


0.95 


1.00 


k B T^ 2) /J 


2.66(1) 


2.55(1) 


2.47(1) 


2.41(1) 


2.36(1) 


2.31(1) 


2.27(1) 


p/vz 


0.019(5) 


0.039(5) 


0.060(4) 


0.094(6) 


0.116(7) 


0.075(4) 


0.057(3) 


z 


1.97(4) 


2.02(3) 


2.10(3) 


2.09(3) 


2.09(6) 


2.20(4) 


2.15(3) 


e 


0.43(3) 


0.21(7) 


0.22(3) 


0.11(4) 


0.16(5) 


0.13(3) 


0.19(4) 


V 


0.07(2) 


0.16(2) 


0.25(2) 


0.39(3) 


0.48(3) 


0.33(2) 


0.25(1) 


r 


0.994455464 


0.998375667 


0.999227272 


0.999226704 


0.99928958 


0.99925661 


0.997206853 



TABLE V. Critical temperature and exponents obtained for different q values for the Metropolis II algorithm. The exponents 
were obtained performing simulations for the estimated critical temperatures. They were based on power laws described in 
Sec. [TTJ The last line we show the r value for the best fits in the second stage (fine scale) 
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t(MCsteps) 



q = 0.70 
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t(MCSteps) 
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FIG. 5. Time correlation of magnetization C(t) 
(M(t)M(0)) for two prescriptions: Metropolis I and II. 



ing from tuq = 1 (ordered state). Similarly, we repeated 
the numerical experiment for N^m = 5 different seeds to 
obtain the uncertainties. In two dimensional systems, the 
slope is (j> — 2/z (see Eq.|4j and so z is calculated accord- 
ing to z = 2/(j) and the uncertainty in z is obtained by 
relation a z = (2/<p 2 )(T ( f,. Here the • denotes the amount 
estimated from Nb in — 5 different seeds. Once z is cal- 
culated, the exponent r\ = 2/3 /V is calculated according 

to r\ = 2(p/vz) ■ z, where (fijvz) was estimated via mag- 
netization decay and z from cumulant F2 . The exponent 
9 was similarly obtained performing N s = 3000 differ- 
ent runs to evolve the time series of correlation C (t) and 
estimating directly the slope in this case. 

Tables HTT1 and IVl show results for the critical exponents 
obtained with the two algorithms. We do not observe 
a monotonic behavior of the critical exponents as func- 
tion of q in either case but on the other hand for both 
cases we cannot assert, for example, that z € [2.09,2.15] 
(Metropolis I) and z e [1.97, 2.20] (Metropolis II) or even 
other exponents do not change for q < 1 which implies 
that we cannot simply extrapolate the critical properties 
from q = 1 to q < 1. 



V. CONCLUSIONS 

In the non-extensive thermostatistics context, we have 
proposed a generalized master equation leading to a gen- 
eralized Metropolis algorithm. This algorithm is local 
and satisfies the detailed energy balance to calculate the 
time evolution of spins systems. We calculate the critical 
temperatures using the generalized Metropolis dynamics, 
via equilibrium and non-equilibrium Monte Carlo simu- 
lations. 

We have obtained the critical parameters performing 
Monte Carlo simulations in two different ways. Firstly, 
we show the phase transitions from curves (M) versus 



UbT / J, considering the magnetization averaging, in equi- 
librium, under different MC steps. Next, we use the 
short time dynamics, via relaxation of magnetization 
from samples initially prepared of ordered or disordered 
states, i.e., time series of magnetization and their mo- 
ments averaged over initial conditions and over different 
runs. 

We have also studied the Metropolis algorithm of 
Refs. [UG3- We show that it does not preserve local- 
ity neither the detailed energy balance in equilibrium. 
While our non-equilibrium simulations corroborate re- 
sults of Refs.fll], EH when we use their extension of the 
Metropolis algorithm (Metropolis I), the exponents and 
critical temperatures obtained are very different when we 
use our prescription (Metropolis II). When the extensive 
case is considered, both methods lead to the same ex- 
pected values. 

Simultaneously, we have developed a methodology to 
refine the determination of the best critical temperature. 
This procedure is based on optimization of the power 
laws of the magnetization function that relaxes from or- 
dered state in log scale, via of maximization of determi- 
nation coefficient of the linear fits. This approach can 
be extended for other spin systems, since their general 
usefulness. 

For a more complete elucidation about existence of 
phase transitions for q ^ 1, we have performed simu- 
lations for small systems MC simulations, recalculating 
the whole lattice energy in each simple spin flip, accord- 
ing to Metropolis I algorithm only to check the varia- 
tions on the critical behavior of the model. Notice that 
this does not apply to Metropolis II algorithm, since it 
has been designed to work as the standard Metropolis 
one. Our numerical results show discontinuities in the 
magnetization, but no finite size scaling, corroborating 
the results of Ref . [HI , which used the broad histogram 
technics to show that no phase transition occurs for q =/= 1 
using Metropolis I algorithm. 

It is important to mention that only Metropolis I 
[Til [r2j shows inconsistence on critical phenomena of 
model since global and local simulation schemes leads to 
different critical properties. Metropolis II overcomes this 
problem since local and global prescriptions are the same 
even for q ^ 1. Broad histogram method works with a 
non-biased random walk that explore the configuration 
space, leading to a phase transition suppression for q =/= 1 
[l3| ■ Nevertheless this algorithm must also be adapted to 
deal with the generalized Boltzmann weight in the same 
way the master equation needed to be modified. This is 
out of the scope of the present paper but this issue will 
be treated in a near future. 
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